#This file creates Figure 1 - Eigenvalues of the correlation matrix and correct 
#classification of the ideal point model in multiple dimensions

rm(list = ls())
library(foreign)
library(wnominate)
library(pscl)
library(readr)

dataset <- read.csv("roll_call_matrix.csv")

rc.matrix <- dataset[,2:ncol(dataset)]
m <- match(rownames(rc.matrix), rownames(dataset))
mat <- rollcall(rc.matrix, yea = 1, nay = 0)

#############################################
#Figure 1a - Dimensionality Analysis#########
#############################################

#Eigenvalues of the correlation matrix
r <- cor(t(rc.matrix) ,use="pairwise")
dim(r)
table(is.na(r))

e <- eigen(r)
lambda <- e$values
n <- dim(e)[1]
plot(1:15, lambda[1:15], ylab = "Eigenvalues", xlab = "Dimension", type = "l", axes = F, main = "Eigenvalues of the Correlation Matrix",
     ylim = c(0, 22), xlim = c(1, 15))
points(1:15,lambda[1:15], pch = 16)
axis(2, at = seq(0, 22, 2), las = 2)
axis(1, at = seq(1, 15, 2))

eigens <- prcomp(na.omit(r), center = T, scale. = T)
prop.var <- summary(eigens)
prop.var <- prop.var$importance[2,]

plot(1:10, prop.var[1:10], type = "l", main = "Dimensionality of Features \n Correlation Matrix", axes = F, ylab = "Proportion of Variance Explained", xlab = "Dimension", ylim = c(0, .8))
points(1:10, prop.var[1:10], pch = 16)
axis(1, at = 1:10, cex.axis = .85)
axis(2, at = seq(0, .8, .1), labels = seq(0, 80, 10), las = 2)

#############################################
#Figure 1b -Correct Predicted Analysis#######
#############################################

#scale that beautiful roll call matrix
set.seed(4227492)
idealpoints <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 1, store.item = T, normalize = T)
y
y
pred.idealpoints <- predict(idealpoints, cutoff=.5)
correct <- pred.idealpoints$overall.percent
#87.18371

#run 2-10d models in ideal to get %correct classified
set.seed(4227492)
idealpoints2 <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 2, store.item = T, normalize = T)
y
y
pred.idealpoints2 <- predict(idealpoints2, cutoff=.5)
correct2 <- pred.idealpoints2$overall.percent

set.seed(4227492)
idealpoints3 <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 3, store.item = T, normalize = T)
y
y
pred.idealpoints3 <- predict(idealpoints3, cutoff=.5)
correct3 <- pred.idealpoints3$overall.percent

set.seed(4227492)
idealpoints4 <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 4, store.item = T, normalize = T)
y
y
pred.idealpoints4 <- predict(idealpoints4, cutoff=.5)
correct4 <- pred.idealpoints4$overall.percent

set.seed(4227492)
idealpoints5 <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 5, store.item = T, normalize = T)
y
y
pred.idealpoints5 <- predict(idealpoints5, cutoff=.5)
correct5 <- pred.idealpoints5$overall.percent

set.seed(4227492)
idealpoints6 <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 6, store.item = T, normalize = T)
y
y
pred.idealpoints6 <- predict(idealpoints6, cutoff=.5)
correct6 <- pred.idealpoints6$overall.percent

set.seed(4227492)
idealpoints7 <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 7, store.item = T, normalize = T)
y
y
pred.idealpoints7 <- predict(idealpoints7, cutoff=.5)
correct7 <- pred.idealpoints7$overall.percent

set.seed(4227492)
idealpoints8 <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 8, store.item = T, normalize = T)
y
y
pred.idealpoints8 <- predict(idealpoints8, cutoff=.5)
correct8 <- pred.idealpoints8$overall.percent

set.seed(4227492)
idealpoints9 <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 9, store.item = T, normalize = T)
y
y
pred.idealpoints9 <- predict(idealpoints9, cutoff=.5)
correct9 <- pred.idealpoints9$overall.percent

set.seed(4227492)
idealpoints10 <- ideal(mat, maxiter = 10000, burnin = 4000, thin = 1, verbose = TRUE, d = 10, store.item = T, normalize = T)
y
y
pred.idealpoints10 <- predict(idealpoints10, cutoff=.5)
correct10 <- pred.idealpoints10$overall.percent

cor.class10 <- c(correct, correct2, correct3, correct4, correct5, 
                 correct6, correct7, correct8, correct9, correct10)
x <- 1:10
plot(x, cor.class10, ylim = c(0,100), pch = 16, 
     xlab = "Dimensions", ylab = "Percent", 
     type = "l", main = "Correct Classification of Votes", axes = F)
points(x,cor.class10, pch = 16)
axis(1, at = 1:10)
axis(2, at = seq(0, 100, 20), las = 2)

